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Abstract 

We propose a new metric to quantify the correlation be- 
tween any two earthquakes. The metric consists of a prod- 
uct involving the time interval and spatial distance be- 
tween two events, as well as the magnitude of the first one. 
According to this metric, events typically are strongly 
correlated to only one or a few preceding ones. Thus a 
classification of events as foreshocks, main shocks or af- 
tershocks emerges automatically without imposing prede- 
fined space-time windows. To construct a network, each 
earthquake receives an incoming link from its most cor- 
related predecessor. The number of aftershocks for any 
event, identified by its outgoing links, is found to be scale 
free with exponent 7 — 2.0(1). The original Omori law 
with p—X emerges as a robust feature of seismicity, hold- 
ing up to years even for aftershock sequences initiated by 
intermediate magnitude events. The measured fat-tailed 
distribution of distances between earthquakes and their 
aftershocks suggests that aftershock collection with fixed 
space windows is not appropriate. 



1 Introduction 

Earthquakes exhibit complex correlations in space, time, 
as well as magnitude |21 El E El El- Sequences of 
earthquakes often appear related to main shocks of large 
magnitude, which are followed in time by nearby smaller 
events. Sometimes, the main shock is also preceded by 
a few intermediate or smaller precursor events. Earth- 
quakes can also cluster as swarms, where the seismic ac- 
tivity is not distinctly associated with a main event. Hu- 
man observation tends toward labeling these events de- 
pending on their relative magnitude and their position 
in the space-time sequence: foreshocks, main shocks and 
aftershocks, respectively. However, in identifying after- 
shocks, it is necessary to distinguish them from what is 
called background seismicity, and to identify their main 
shock. Although an observation by eye of the evolving 
seismic situation can support a classification, a precise 
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Figure 1: Schematic examples of space-time windows 
used to collect aftershocks: the usual rectangular or con- 
vex window (dashed line) and our hyperbolic, concave 
window (shaded region). 

label for each event may be intrinsically impossible. 

In the most popular approach, aftershocks are collected 
by counting all events within a fixed space-time win- 
dow 7 , 8 , 9i ilOj following a pre-assigned main event (see 
Fig. However, this method does not define the prob- 
ability that an event thereby collected is actually corre- 
lated to the main event under consideration. Maybe more 
importantly, one does not know whether the predefined 
space-time windows are too large or too small for min- 
imizing errors in the procedure. A more subtle issue is 
to define aftershocks of aftershocks. If an aftershock can 
have more than one preceding large event, which of these 
should be regarded as the most important or correlated 
one? 

A quantitative metric of the correlation between any 
two earthquakes, or the extent to which one can be consid- 
ered an aftershock of another, may be crucial for solving 
these problems, and for developing a better understand- 
ing of seismicity. A reliable metric should include known 
statistical properties. One is the Gutenberg-Richter (G- 
R) distribution 4 for the number of earthquakes of mag- 
nitude TO in a seismic region. 
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with b usually « 1. Another is the fractal appearance of 
earthquake epicenters PElEji with fractal dimension df. 
Thus, the average number of earthquakes of magnitude 
within an interval Am of to, occurring in an area of radius 
r over a time interval r, is 
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where C is a constant depending on the overall seismicity 
in the region under consideration. 

For any earthquake j in the seismic region, looking 
backward in time, how many earthquakes of magnitude 
within an interval Am of m would be expected to have 
occurred within a time interval t, and within a distance /, 
of that specific event? In fact, an n value can be defined 
between any two events i and j occurring in the sequence 
at times Ti and Tj, with Ti < Tj. If we take the mag- 
nitude rui of the z-th event, the spatial distance I = hj 
between the two earthquake epicenters, and the time in- 
terval t = tij = Tj — Ti, the expected number of events of 
magnitude within Am of m^ occurring in the particular 
space-time domain bounded by events j and i is 

= Ctl'^f Am 10-''"'^ . (3) 

Note that the domain appearing in Eq. O is selected by 
the particular history of the seismic activity in the region 
and not predefined by any observer. 

Of all the earthquakes preceding j, the most unlikely 
to occur according to Eq.Elis earthquake i* such that riij 
is minimized when i = i*. However, earthquake i* actu- 
ally occurred relative to j, even though it was the least 
likely to have done so. Therefore, i* must be the event 
to which earthquake j is most correlated. In general, if 
riij is extremely small, then the correlation between j 
and i is very strong, and vice versa. By this argument, 
the degree of correlation between any two earthquakes 
i and j is inversely proportional to Uij . Since the space- 
time-magnitude scales are selected by the actual sequence 
of events, the variables riy can be considered to be self- 
organizing tags of the underlying physical process gov- 
erning seismicity. Note that singularities are eliminated 
by taking a small scale cutoff in time (here t^ain = 180 
sec) and a minimum spatial resolution (here ^juin = 100 
meters) . 

The metric defined by Eq. O allows various classifica- 
tions of aftershocks. Therefore, the question of which 
is the better candidate to be the foreshock of an event 
can be quantitatively decided. Hierarchical clusters of 
earthquakes emerge, in which the biggest event in the 
cluster is called the main event, but where possibly later 
aftershocks create their own sequences of aftershocks, 
whenever they are able to "steal" aftershocks from the 
main event, and so on for further generations of after- 
shocks. Nevertheless, earthquakes are automatically col- 
lected into hierarchically self-organized clusters, without 
any special pre-analysis of single event properties. 

In the language of modern complex network theory 
|12[ I13| , what we achieve is a time-oriented growing net- 
work where nodes (earthquakes) have internal variables 
(magnitude, occurrence time, and location), and links be- 
tween the nodes carry a weight (the metric Uij) and are 
directed according to the time orientation, from the older 
to the newer nodes. Empirically, we find that both the 
distribution of outgoing links and the cluster size distri- 
bution are scale free. Due to the continuous nature of the 



link variable, riij, no event is a priori purely an aftershock 
or a main shock. However, due to the broad distribution 
of Uij observed, main shocks and aftershocks emerge as 
extreme limits of a continuous spectrum of the extent to 
which any given event can be considered to be a precursor 
or aftershock of other events in the sequence. 

Our approach was inspired by a recent analysis of earth- 
quake waiting times by Bak et al. [HI Ej- They intro- 
duced a space-time-magnitude scaling variable that al- 
lows a data collapse of the distribution of waiting times 
between subsequent earthquakes larger than a specified 
magnitude, occurring within grid cells of a specified size, 
covering non-overlapping areas of the Earth. Also, Abe 
and Suzuki found scale free networks for earthquakes in 
a completely different context, where nodes representing 
these grid cells were linked when subsequent earthquakes 
occurred in them |15| . However, neither of these works 
quantified the correlation between an arbitrary pair of 
earthquakes, or dealt with the subject of aftershock iden- 
tification. 

2 Data and parameters 

The catalog we have analyzed is maintained by 
the Southern California Earthquake Data Center 
(it can be downloaded from the SCEDC web site 
http: / /www. scecdc.scec.org/ ftp/catalogs / scsn) , for 
which Am = 0.1. It is considered to be complete 
for events with m > 2. We use data ranging from 
January 1, 1984 to December 31, 2000. In order to work 
with a well defined ensemble, a lower threshold on the 
magnitude is introduced: events with magnitude smaller 
than m< are discarded. For each event, its position 
i in the sequence is used as a label, and we record 
the magnitude rrii, the occurrence time Ti (measured 
in seconds from midnight of the first day), and the 
latitude and longitude of the epicenter (converted to 
angles measured in radians, 0i and (j)i respectively). 
The distance between two events i and j is then mea- 
sured as the arc length on the Earth's surface, kj = 
_Ro arccos[sin(0i) sin(6'j) + cos(0i) cos(0j) cos(0i — (f)j)], 
where the Earth radius is Rq = 6.3673 x 10^ meters. 

The 6-value of the G-R law is 6 ~ 0.95 for this data 
set, while df ~ 1.6 was found by Corral lA using a box 
counting procedure. It is consistent with the correlation 
dimension we measure for most of our clusters. However, 
many of the statistical results we find are not sensitive to 
the precise value oi df or b. 

With these units and values, the constant C can be 
estimated using Eq. El However, a precise evaluation of 
C is not possible, because n is the mean of a variable with 
huge variations in space and time. We have measured n 
for several circular windows well inside the zone covered 
by the catalog, finding C < 10~^. For simplicity, our 
choice in this paper is C = 10^^. Most of our results are 
insensitive to the precise value of C because we focus on 
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relative, rather than absolute correlations between a pair 
of events. Throughout this paper we use, unless otherwise 
stated, the above mentioned values, and a lower threshold 
TO< = 2.5. 

To simplify notation, we denote the probability distri- 
bution of a generic quantity q as P{q). On finding distri- 
butions decaying as power laws, a clearer result appears 
by binning the values of P{q) in properly normalized bins 
of a width that grows geometrically with q. 

3 Method 

In the simplest implementation of the network, each new 
earthquake j is attached with a single link to the previ- 
ous earthquake in the sequence that minimizes n^j, with 
a weight denoted as n*. Hence, each link carries the ex- 
tremal n* for the added node j relative to all previous 
nodes, and globally one obtains a growing directed tree. 
Links with small n* indicate a stronger correlation be- 
tween the emitting node and the receiving one, and are 
expected to identify events normally classified as after- 
shocks. Weak links with large n* arise when none of the 
previous events are sufficiently strong, and close in space 
and time to event j. Thus, the strength of the link to an 
event j is inversely proportional to n* . 

A natural decomposition of the network into clusters is 
achieved by then removing all weak links where n* > ric, 
and Tic is a link threshold value. The correlated events 
are reliably detected when ric is less than one but not 
extremely small. In the latter case, correlated events de- 
tach, and a very fragmented network appears. For large 
ric some uncorrelated events make links, and a giant clus- 
ter appears. The resulting space-time windows are con- 
cave (see Figure^ and Conclusions), at variance with the 
convex windows usually used. 

4 Results 

A part of the network constructed using this method 
is shown in Fig. El Hierarchically organized clusters of 
earthquakes emerge, where the links join aftershocks with 
their most correlated predecessor. 

In order to quantitatively assess the properties of this 
network, we start by analyzing the distribution of link 
weights P{n*). This distribution exhibits power law be- 
havior with a slope ~ —1 up to a cutoff, as shown in Fig. 13 
Thus, the distribution of link strengths, S — l/n*, is also 
a power law, P{S) ~ 1/S. Such a broad, continuous 
distribution, without particular characteristic peaks, in- 
dicates that a division of earthquakes into rigid classes is 
intrinsically impossible. Instead, a continuum of possibil- 
ities ranges from clear aftershocks, which have an incom- 
ing link with small n*, to events that are independent, 
with an incoming link of large n*, but may emit many 
outgoing links with small n*, and would be called main 
shocks. 




Figure 2: Scale free earthquake network around Landers 
epicenter (red cluster) and Hector Mine epicenter (blue 
cluster) . Colors fade with the aftershock generation, from 
darker to lighter within each cluster. Note that the big 
event following the Landers earthquake, giving rise to its 
own (orange) cluster of aftershocks at (—116.50', 34.10'), 
is not a first generation aftershock, since it has no link 
from Landers. Here m< — 4, and ric = 10^^. 

The resulting network of earthquakes is scale free. The 
number of aftershocks of an earthquake is equal to the 
number, fc, of outgoing links from the node represent- 
ing that event. In the language of network theory, this 
is called the out-degree of the node. Fig. 0] shows that 
earthquakes in Southern California form a scale free net- 
work, with an out-degree distribution scaling over more 
than three decades, with an index 7 — 2.0(1). 

Recently, many scale free networks with P{k) ~ k'"', 
have been discovered ^1 El in a broad variety of con- 
texts. These include the Internet ^HI, the world-wide 
web |17| . protein interaction and genetic regulatory net- 
works ^1 , and the solar coronal magnetic field . 
Both the coronal field and the earthquake network share 
the property that the strength of a link between a pair of 
nodes is also scale free. 

Lowering the link threshold ric from infinity, the fully 
connected network breaks into clusters, in a percolation- 
like transition from a giant cluster to a finite cluster 
regime. We estimate the transition to take place between 
link thresholds ric = 10^^ and Uc = 10^^. This esti- 
mate is obtained by examining the distribution of cluster 
sizes, N, which is the total number of earthquakes in a 
cluster, as a function of ric (see Fig. jSJ. Near the transi- 
tion, the cluster size distribution also appears to be scale 
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Figure 3: The distribution of link weights, n*, for se- 
quences of different temporal duration. An average over 
all non-overlapping time intervals of the same duration 
is shown. The power law behavior is stable to variations 
in the duration. However, the cutoff moves to smaller n* 
on increasing the measurement time interval as weakly 
linked earthquakes find more correlated predecessors fur- 
ther in the past. The vertical dotted line represents the 
estimated transition point, n^, for the giant cluster. The 
straight dashed line has a slope -1. 
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Figure 4: The degree distribution of the network of earth- 
quakes and aftershocks. The out-degree, k, is the number 
of aftershocks linked to an earthquake. The introduction 
of a threshold ric does not alter the behavior. The dashed 
line has slope —2, indicating a scale free degree distribu- 
tion, P{k) - fc~T with 7^2. 

free, P{N) ~ A^~^'^*^^^ Furthermore, a scaling regime 
exists for a wide range of link thresholds, indicating a rel- 
ative insensitivity to a sharp separation between what are 
considered to be correlated and uncorrelated events. For 
clarity, we use the value = 10~^ to locate the tran- 
sition point where the giant cluster emerges. This value 
is consistent with our ansatz, Eq. |31 which requires that 
correlated events have n values significantly less than one. 
Network clusters constructed with therefore only link 
strongly correlated events. 
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Figure 5: Cluster size distribution for different link 
thresholds. At large ric, a giant cluster exists that is 
well separated in size from some remaining small ones. 
Between ric — 10^^ and ric = 10~^, an apparently contin- 
uous transition occurs where the finite cluster distribution 
extends out toward the giant cluster, and the distribution 
of cluster sizes exhibits power law behavior. The straight 
line has a slope -1.7. Symbols are (*ric = 10^; Oric = 
10; On, = 1; +n, ^ lO^i; Q = 10^')- 



In Fig. 13 we study the effect of changing the temporal 
span of the catalog on the distribution of link weights. 
The power law behavior for strong links is stable, cer- 
tainly up to ric- However, the cutoff in P{n*) for weak 
links decreases to smaller n* values, when earthquakes 
can link to events at further distance in the past. For 
an ideal "infinite" catalog, we conjecture that the cutoff 
value cannot be less than . Indeed, below the transition 
point a finite fraction of events stop having any correlated 
predecessors. 

We define the link length, /, as the distance between 
the epicenter of an aftershock and its linked predecessor. 
The distribution of link lengths depends on the magni- 
tude m of the predecessor, being on average greater for 
larger m. Dividing the link length distribution into classes 
depending on the magnitude of the predecessor, Pjn{l), a 
maximum in the distribution occurs, which shifts to larger 
I on increasing m, as shown in Fig. |3 This behavior is 
consistent with using larger space-time windows to collect 
aftershocks from larger events. 

However, the distribution of link lengths exhibits no 
cutoff at large distances, but rather decays slowly as a 
power law with /, up to the linear extent of the seismic 
region covered by the catalog. The different distributions 
are consistent with a scaling ansatz: 



~ 10-''™F(Z/10'"'") 



(4) 



where I is measured in kilometers, a ~ 0.4, and F{x) is 
a scaling function. The tail of the scaling function is a 
power law; i.e. F{x) ^ with A w 2 for a; ^ 1. A data 
collapse using this ansatz is shown in the inset to Fig. 
Such a slow decay at large distances calls into question 
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Figure 6: Link length distribution for different mag- 
nitudes of the emitting earthquake, at n^. The length 
at maximum grows with magnitude roughly as /max ~ 
10°'^™, but the distributions have a fat tail, extending 
up to hundreds of kilometers even for intermediate mag- 
nitude events. These distributions are consistent with 
an hierarchical organization of events, where big earth- 
quakes preferentially link at long distance with interme- 
diate ones, which in turn link to more localized after- 
shocks, and so on. Inset: Distributions rescaled according 
to Eq.Hwith a = 0.4. 

the use of sharply defined space windows for collecting 
aftershocks, as already pointed out by Ogata [21] ■ 

Figure H shows the rate of aftershocks for the Landers, 
Hector Mine, and Northridge events. Aftershocks occur- 
ring at time t after one of these events are binned into 
geometrically increasing time intervals. The number of 
aftershocks in each bin is then divided by the temporal 
width of the bin to obtain a rate of earthquakes per sec- 
ond. The same procedure is applied to each remaining 
event, not aftershocks of these three. An average is made 
for the rate of aftershocks linked to events having a mag- 
nitude within an interval Am of to. Figure also shows 
the averaged results for m = 3 (1710 events), m — A (161 
events), m = 5 (28 events) and to = 5.9 (4 events). 

The collection of aftershocks linked to earthquakes of 
all magnitudes is one of the main results of our method. 
Even intermediate magnitude events can have aftershocks 
that persist up to years. Earthquakes of all magnitudes 
have aftershocks which decay according to the Omori 
law 0122, 

K 

'^(*) ^ — r~r ' fo^' ^ < ^cutoff (5) 

c + t 

where c and K are constant in time, but depend on the 
magnitude to |22[I23) of the earthquake. We find that the 
Omori law persists up to a time ^cutoff that also depends 
on TO as well as the link threshold, Uc- Estimates of the 
cutoff times for are tcutoff ~ 3 months for m = 3, and 
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Figure 7: The Omori law for aftershock rates. These 
rates are measured for aftershocks linked to earthquakes 
of different magnitudes, to, using n^. For each magnitude, 
the rate is consistent with the original Omori law, Eq. 
up to a cutoff time that depends on to. As guides to the 
eye, dashed lines represent a decay ~ 1/t. 

i^^jQg « 1 year for to = 4. For larger magnitudes, it is 
difficult to distinguish icutoff from the temporal duration 
of the data set. 

The Omori law for aftershocks emerges as a result of 
our analysis, although it is not part of the original ansatz, 
Eq. 13 used to define aftershocks. It has been extensively 
investigated over decades, together with its modified ver- 
sion 122] involving a scaling ^ t^P. The data shown 
in Figure [3 is consistent with the original Omori result, 
p — 1, for aftershocks of earthquakes of all magnitudes, 
once second and further generations of aftershocks are 
excluded. Our result is also consistent with theoretical 
studies on stick-slip motion ,24, ,25,, which suggest p « 1. 

5 Discussion 

Convex space time windows have been used since the 
1970's iJ, |8i i9| ilOj , often with the size of the window de- 
termined by the main shock magnitude. The performance 
of this procedure is satisfactory for large earthquakes, al- 
though fixed window sizes may omit relevant aftershocks. 
Nevertheless, as a shortcoming, it can lead to distortions 
if many large aftershocks occur. In this case, nothing can 
be said on the "ownership" of further aftershocks. 

Different approaches to the problem of aftershocks col- 
lection were proposed by several authors, sometimes with 
the aim to cure the former shortcomings. For a review 
see |2E1- Our method has some similarities with these 
approaches. For instance, Frohlich and Davis collected 
earthquakes in clusters j27j by means of a different linking 
procedure. However, their analysis was done using a met- 
ric of the form ^ y/P + const t^, which does not take into 
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account the magnitude of events, and has a space-time 
form at variance with measured earthquake correlations. 

Maximum hkehhood methods [211123 1 in the context of 
seismicity, usually start with an ansatz on the law govern- 
ing aftershocks, typically the modified Omori law. It is 
further assumed that seismicity is a non-stationary Pois- 
son process. Using a likelihood analysis with space, time 
and magnitude, Ogata compared several forms of after- 
shocks distance distributions |2H . and showed that an 
aftershock rate of the form 

was the most appropriate among his choices (ct, a, p and 
fi are constant, while ci{'m) is scaling with the magnitude 
TO of the main shock). Hence, he also concluded that 
fixed space windows were not the best choice. Indeed, 
our metric variable n in Eq. |31 somewhat resembles his 
form of V. 

However, our method is simpler to implement than 
likelihood methods. Furthermore, it does not require an 
ansatz on the validity of the modified Omori law, or on 
the type of statistical process that describes seismicity. 
Instead, the original Omori law is found as a result of 
our analysis. In addition, the physical argument leading 
to the variable also fixes the parameters in its defini- 
tion, without the need to evaluate them by maximizing a 
likelihood. 

One could object that the values of 6 and/or df can de- 
pend on the region of the Earth being considered, or may 
fluctuate depending on the specific fault zone being stud- 
ied. However, the statistical results we find, as shown in 
Figs. 3-7, are remarkably robust to variations in either of 
these parameters. Varying df over a wide range, from 1 to 
3 (using df > 2 requires the introduction of event depths, 
see below) does not alter considerably the distribution of 
outgoing links, which retains its power law behavior with 
index 7 « 2. The distribution of link weights, n*, is even 
more insensitive to variations of b and d f . Also the Omori 
law with p « 1, shown in Fig. 7, does not depend sensi- 
bly on the parameters, and holds for aftershocks linked 
to earthquakes of all magnitudes. 

The crust of the Earth has a finite width 20 Km 
in California) in which events take place according to a 
"three-dimensional" fractal distribution, involving their 
depth. It is believed that there is a qualitative difference 
between small earthquakes and large ones, the former pro- 
ducing ruptures smaller that the crust width T. Hence, 
our arguments may need to be corrected at distances of 
the order of tens of kilometers. We have computed spatial 
distances through the three-dimensional Euclidean met- 
ric distance, using an appropriately revised df in Eq. 13 
No significant departures from the results leading to our 
present conclusions were found. 

The introduction of more than one correlated predeces- 
sor for an event will be the subject of a future investiga- 
tion. This is accomplished by attaching links between all 



earthquake pairs where riy < ?ic- In this case, a general 
network, which is not tree-like emerges. The clustering 
of earthquakes could then be quantified in terms of the 
clustering coefficient of the nodes in the network [301 E| ■ 
In our view, an earthquake network with nodes having 
multiple incoming links represents a second order model- 
ing of seismicity, the first being the simple tree structure 
we have presented here. 

6 Conclusions 

We have introduced a metric to determine correlations 
between earthquakes that takes into account known sta- 
tistical properties of seismicity. By means of an appeal- 
ingly simple yet quantifiable procedure, networks of earth- 
quakes and aftershocks emerge, where the number of af- 
tershocks linked to any event is scale free with an index 
7^2. The metric is constructed by looking backward 
in time from any particular event and calculating an ex- 
pected number of events that would occur, compared to 
events that actually occurred. If this ratio is significantly 
less than one, then the preceding event is correlated with 
the particular one. This is reminiscent of Kierkegaard's 
adage that life must be lived forward, but can only be 
understood backward. 

Due to the form of the metric n measuring correla- 
tions, larger earthquakes collect aftershocks from larger 
space-time windows. From Eq. |31 these windows have 
a spatial radius varying with time as ri(T) = [ndT — 
T,)-'^10^'""]^/'^f . They span an hyperbolic space-time re- 
gion (see Fig.^, which is at variance with the usual "rect- 
angular" or convex windows, of constant radius up to a 
finite time. In our method, at early times after an earth- 
quake, its aftershock collection window is wider in space 
than it is at later times. 

To our knowledge, the idea that an earthquake can be 
correlated to an event very far away, if it occurs shortly 
after it, is new and certainly unconventional. But it is 
consistent with the hypothesis that seismicity is a self- 
organized critical phenomenon |31l I32L In that case, 
some locations may be "on the edge of giving an earth- 
quake" (or toppling, according to the sandpile paradigm), 
and even a small perturbation from an event far away 
could trigger them. However, we do not necessarily as- 
cribe the correlations measured here to represent a usual 
cause and effect relationship. In the sandpile paradigm 
a completely insignificant event, like adding one grain of 
sand to an enormous pile, can trigger an arbitrarily large 
avalanche involving the whole system. Indeed, seismic- 
ity as one hierarchically correlated self-organized critical 
process, generates the scale free network of earthquakes 
and aftershocks. 

Our results also suggest that modern network theory 
may be a useful and illuminating way to approach the 
complexities of seismicity, including perhaps problems re- 
lated to prediction. Our metric and network construction 
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may also have applications to other phenomena with in- 
termittent bursts such as, for instance, solar flares or even 
turbulence. 
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